/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
	This file is part of foam-extend.

	foam-extend is free software: you can redistribute it and/or modify it
	under the terms of the GNU General Public License as published by the
	Free Software Foundation, either version 3 of the License, or (at your
	option) any later version.

	foam-extend is distributed in the hope that it will be useful, but
	WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
	General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "cyclicPolyPatch.H"
#include "addToRunTimeSelectionTable.H"
#include "polyBoundaryMesh.H"
#include "polyMesh.H"
#include "demandDrivenData.H"
#include "OFstream.H"
#include "patchZones.H"
#include "matchPoints.H"
#include "EdgeMap.H"
#include "foamTime.H"
#include "RodriguesRotation.H"
#include "standAlonePatch.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
	defineTypeNameAndDebug(cyclicPolyPatch, 0);

	addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, word);
	addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, dictionary);


template<>
const char* NamedEnum<cyclicPolyPatch::transformType, 3>::names[] =
{
	"unknown",
	"rotational",
	"translational"
};

const NamedEnum<cyclicPolyPatch::transformType, 3>
	cyclicPolyPatch::transformTypeNames;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Foam::label Foam::cyclicPolyPatch::findMaxArea
(
	const pointField& points,
	const faceList& faces
)
{
	label maxI = -1;
	scalar maxAreaSqr = -GREAT;

	forAll(faces, faceI)
	{
		scalar areaSqr = magSqr(faces[faceI].normal(points));

		if (areaSqr > maxAreaSqr)
		{
			maxAreaSqr = areaSqr;
			maxI = faceI;
		}
	}
	return maxI;
}


void Foam::cyclicPolyPatch::calcTransforms() const
{
	if (size())
	{
		const pointField& points = this->points();

		// Determine geometric quantities on the two halves
		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		primitivePatch half0
		(
			SubList<face>
			(
				*this,
				size()/2
			),
			points
		);

		pointField half0Ctrs(calcFaceCentres(half0, half0.points()));
		scalarField half0Tols
		(
			calcFaceTol(half0, half0.points(), half0Ctrs)
		);

		primitivePatch half1
		(
			SubList<face>
			(
				*this,
				size()/2,
				size()/2
			),
			points
		);
		pointField half1Ctrs(calcFaceCentres(half1, half1.points()));

		// Dump halves
		if (debug)
		{
			fileName fvPath(boundaryMesh().mesh().time().path()/"VTK");

			mkDir(fvPath);
			fileName nm0(fvPath/name() + "_half0_faces");
			Pout<< "cyclicPolyPatch::calcTransforms : Writing half0"
				<< " faces to file " << nm0 << endl;
			half0.writeVTK(nm0, half0, points);

			fileName nm1(fvPath/name() + "_half1_faces");
			Pout<< "cyclicPolyPatch::calcTransforms : Writing half1"
				<< " faces to file " << nm1 << endl;
			half1.writeVTK(nm1, half1, points);
		}

		vectorField half0Normals(half0.size());
		vectorField half1Normals(half1.size());

		scalar maxMatchError = 0;
		label errorFace = -1;

		for (label facei = 0; facei < size()/2; facei++)
		{
			half0Normals[facei] = operator[](facei).normal(points);
			label nbrFacei = facei + size()/2;
			half1Normals[facei] = operator[](nbrFacei).normal(points);

			scalar magSf = mag(half0Normals[facei]);
			scalar nbrMagSf = mag(half1Normals[facei]);
			scalar avSf = (magSf + nbrMagSf)/2.0;

			if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL)
			{
				// Undetermined normal. Use dummy normal to force separation
				// check. (note use of sqrt(VSMALL) since that is how mag
				// scales)
				half0Normals[facei] = point(1, 0, 0);
				half1Normals[facei] = half0Normals[facei];
			}
			else if
			(
				mag(magSf - nbrMagSf)/avSf
			 > polyPatch::matchTol_()
			)
			{
				// Error in area matching.  Find largest error
				maxMatchError =
					Foam::max(maxMatchError, mag(magSf - nbrMagSf)/avSf);
				errorFace = facei;
			}
			else
			{
				half0Normals[facei] /= magSf;
				half1Normals[facei] /= nbrMagSf;
			}
		}


		// See if transformation is prescribed
		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		if (transform_ == ROTATIONAL)
		{
			if (debug)
			{
				Info<< "Prescribed transform: "
					<< transformTypeNames[ROTATIONAL]
					<< ".  Calculating transforms using Rodrigues Rotation.  "
					<< "Axis = " << rotationAxis_
					<< " centre = " << rotationCentre_
					<< " angle = " << rotationAngle_ << endl;
			}

			forwardT_ = tensorField
			(
				1,
				RodriguesRotation(rotationAxis_,  rotationAngle_)
			);

			reverseT_ = tensorField
			(
				1,
				RodriguesRotation(rotationAxis_, -rotationAngle_)
			);
		}
		else if (transform_ == TRANSLATIONAL)
		{
			if (debug)
			{
				Info<< "Prescribed transform: "
					<< transformTypeNames[TRANSLATIONAL] << endl;

				forwardT_.setSize(0);
				reverseT_.setSize(0);
			}
		}
		else
		{
			if (debug)
			{
				Info<< "Unknown transform: calculating transform from "
					<< "geometry" << endl;
			}

			// Calculate transformation tensors
			calcTransformTensors
			(
				half0Ctrs,
				half1Ctrs,
				half0Normals,
				half1Normals,
				half0Tols
			);

			// Check transformation tensors
			if (forwardT_.size() > 1 || reverseT_.size() > 1)
			{
				SeriousErrorIn
				(
					"void cyclicPolyPatch::calcTransforms() const"
				)   << "Transformation tensor is not constant for the cyclic "
					<< "patch " << name()
					<< ".  Please reconsider your setup and definition of "
					<< "cyclic boundaries." << endl;
			}
		}

		// Dump transformed first half
		if (debug)
		{
			if (reverseT_.size() > 0)
			{
				fileName fvPath(boundaryMesh().mesh().time().path()/"VTK");

				pointField transformPoints = half0.localPoints();

				forAll (transformPoints, pointI)
				{
					transformPoints[pointI] =
					    Foam::transform(reverseT_[0], transformPoints[pointI]);
				}

				standAlonePatch transformHalf0
				(
					half0.localFaces(),
					transformPoints
				);

				fileName nm2(fvPath/name() + "_transform_half0_faces");
				Pout<< "cyclicPolyPatch::calcTransforms : Writing "
					<< "transform_half0 faces to file " << nm2 << endl;

				transformHalf0.writeVTK(nm2, transformHalf0, transformPoints);
			}
		}

		// Check for error in face matching
		if (maxMatchError > polyPatch::matchTol_())
		{
			label nbrFacei = errorFace + size()/2;
			scalar magSf = mag(half0Normals[errorFace]);
			scalar nbrMagSf = mag(half1Normals[errorFace]);
			scalar avSf = (magSf + nbrMagSf)/2.0;

			FatalErrorIn
			(
				"cyclicPolyPatch::calcTransforms()"
			)   << "face " << errorFace
				<< " area does not match neighbour "
				<< nbrFacei << " by "
				<< 100*mag(magSf - nbrMagSf)/avSf
				<< "% -- possible face ordering problem." << endl
				<< "patch:" << name()
				<< " my area:" << magSf
				<< " neighbour area:" << nbrMagSf
				<< " matching tolerance:" << polyPatch::matchTol_()
				<< endl
				<< "Mesh face:" << start() + errorFace
				<< " vertices:"
				<< IndirectList<point>(points, operator[](errorFace))()
				<< endl
				<< "Neighbour face:" << start() + nbrFacei
				<< " vertices:"
				<< IndirectList<point>(points, operator[](nbrFacei))()
				<< endl
				<< "Other errors also exist, only the largest is reported. "
				<< "Please rerun with cyclic debug flag set"
				<< " for more information."
				<< abort(FatalError);
		}

		if (debug)
		{
			// Check the transformation
			scalar maxDistance = 0;
			scalar maxRelDistance = 0;

			forAll (half0Ctrs, faceI)
			{
				if (reverseT_.size())
				{
					maxDistance =
					    Foam::max
					    (
					        maxDistance,
					        mag
					        (
					            Foam::transform(reverseT_[0], half0Ctrs[faceI])
					          - half1Ctrs[faceI]
					        )
					    );

					maxRelDistance =
					    Foam::max
					    (
					        maxRelDistance,
					        mag
					        (
					            Foam::transform(reverseT_[0], half0Ctrs[faceI])
					          - half1Ctrs[faceI]
					        )
					       /(
					           mag(half1Ctrs[faceI] - half0Ctrs[faceI])
					         + SMALL
					        )
					    );
				}
				else
				{
					// Disable checking for translational distance
					// HJ, 13/Jan/2011
//                     maxDistance =
//                         Foam::max
//                         (
//                             maxDistance,
//                             mag
//                             (
//                                 half0Ctrs[faceI]
//                               - half1Ctrs[faceI]
//                             )
//                         );

//                     maxRelDistance =
//                         Foam::max
//                         (
//                             maxRelDistance,
//                             mag
//                             (
//                                 half0Ctrs[faceI]
//                               - half1Ctrs[faceI]
//                             )
//                            /(
//                                mag(half1Ctrs[faceI] - half0Ctrs[faceI])
//                              + SMALL
//                             )
//                         );
				}
			}

			// Check max distance between face centre and
			// transformed face centre
			if (maxRelDistance > sqrt(polyPatch::matchTol_()))
			{
				SeriousErrorIn
				(
					"void cyclicPolyPatch::calcTransforms()"
				)   << "Relative difference in transformation for patch "
					<< name() << " = " << 100*maxRelDistance
					<< "%, abs distance = " << maxDistance
					<< ".  Please check the definition of the transformation"
					<< endl;
			}
		}
	}
}


// Get geometric zones of patch by looking at normals.
// Method 1: any edge with sharpish angle is edge between two halves.
//           (this will handle e.g. wedge geometries).
//           Also two fully disconnected regions will be handled this way.
// Method 2: sort faces into two halves based on face normal.
bool Foam::cyclicPolyPatch::getGeometricHalves
(
	const primitivePatch& pp,
	labelList& half0ToPatch,
	labelList& half1ToPatch
) const
{
	// Calculate normals
	const vectorField& faceNormals = pp.faceNormals();

	// Find edges with sharp angles.
	boolList regionEdge(pp.nEdges(), false);

	const labelListList& edgeFaces = pp.edgeFaces();

	label nRegionEdges = 0;

	forAll(edgeFaces, edgeI)
	{
		const labelList& eFaces = edgeFaces[edgeI];

		// Check manifold edges for sharp angle.
		// (Non-manifold already handled by patchZones)
		if (eFaces.size() == 2)
		{
			if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
			{
				regionEdge[edgeI] = true;

				nRegionEdges++;
			}
		}
	}


	// For every face determine zone it is connected to (without crossing
	// any regionEdge)
	patchZones ppZones(pp, regionEdge);

	if (debug)
	{
		Pout<< "cyclicPolyPatch::getGeometricHalves : "
			<< "Found " << nRegionEdges << " edges on patch " << name()
			<< " where the cos of the angle between two connected faces"
			<< " was less than " << featureCos_ << nl
			<< "Patch divided by these and by single sides edges into "
			<< ppZones.nZones() << " parts." << endl;


		// Dumping zones to obj files.

		labelList nZoneFaces(ppZones.nZones());

		for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
		{
			OFstream stream
			(
				boundaryMesh().mesh().time().path()
			   /name()+"_zone_"+Foam::name(zoneI)+".obj"
			);
			Pout<< "cyclicPolyPatch::getGeometricHalves : Writing zone "
				<< zoneI << " face centres to OBJ file " << stream.name()
				<< endl;

			labelList zoneFaces(findIndices(ppZones, zoneI));

			forAll(zoneFaces, i)
			{
				writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
			}

			nZoneFaces[zoneI] = zoneFaces.size();
		}
	}


	if (ppZones.nZones() == 2)
	{
		half0ToPatch = findIndices(ppZones, 0);
		half1ToPatch = findIndices(ppZones, 1);
	}
	else
	{
		if (debug)
		{
			Pout<< "cyclicPolyPatch::getGeometricHalves :"
				<< " falling back to face-normal comparison" << endl;
		}
		label n0Faces = 0;
		half0ToPatch.setSize(pp.size());

		label n1Faces = 0;
		half1ToPatch.setSize(pp.size());

		// Compare to face 0 normal.
		forAll(faceNormals, faceI)
		{
			if ((faceNormals[faceI] & faceNormals[0]) > 0)
			{
				half0ToPatch[n0Faces++] = faceI;
			}
			else
			{
				half1ToPatch[n1Faces++] = faceI;
			}
		}
		half0ToPatch.setSize(n0Faces);
		half1ToPatch.setSize(n1Faces);

		if (debug)
		{
			Pout<< "cyclicPolyPatch::getGeometricHalves :"
				<< " Number of faces per zone:("
				<< n0Faces << ' ' << n1Faces << ')' << endl;
		}
	}

	if (half0ToPatch.size() != half1ToPatch.size())
	{
		fileName casePath(boundaryMesh().mesh().time().path());

		// Dump halves
		{
			fileName nm0(casePath/name()+"_half0_faces.obj");
			Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
				<< " faces to OBJ file " << nm0 << endl;
			writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());

			fileName nm1(casePath/name()+"_half1_faces.obj");
			Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
				<< " faces to OBJ file " << nm1 << endl;
			writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
		}

		// Dump face centres
		{
			OFstream str0(casePath/name()+"_half0.obj");
			Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
				<< " face centres to OBJ file " << str0.name() << endl;

			forAll(half0ToPatch, i)
			{
				writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
			}

			OFstream str1(casePath/name()+"_half1.obj");
			Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
				<< " face centres to OBJ file " << str1.name() << endl;
			forAll(half1ToPatch, i)
			{
				writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
			}
		}

		SeriousErrorIn
		(
			"cyclicPolyPatch::getGeometricHalves"
			"(const primitivePatch&, labelList&, labelList&) const"
		)   << "Patch " << name() << " gets decomposed in two zones of"
			<< "inequal size: " << half0ToPatch.size()
			<< " and " << half1ToPatch.size() << endl
			<< "This means that the patch is either not two separate regions"
			<< " or one region where the angle between the different regions"
			<< " is not sufficiently sharp." << endl
			<< "Please adapt the featureCos setting." << endl
			<< "Continuing with incorrect face ordering from now on!" << endl;

		return false;
	}
	else
	{
		return true;
	}
}


// Given a split of faces into left and right half calculate the centres
// and anchor points. Transform the left points so they align with the
// right ones.
void Foam::cyclicPolyPatch::getCentresAndAnchors
(
	const primitivePatch& pp,
	const faceList& half0Faces,
	const faceList& half1Faces,
	pointField& ppPoints,
	pointField& half0Ctrs,
	pointField& half1Ctrs,
	pointField& anchors0,
	scalarField& tols
) const
{
	// Get geometric data on both halves.
	half0Ctrs = calcFaceCentres(half0Faces, pp.points());
	anchors0 = getAnchorPoints(half0Faces, pp.points());
	half1Ctrs = calcFaceCentres(half1Faces, pp.points());

	vector n0 = vector::zero;
	vector n1 = vector::zero;
	switch (transform_)
	{
		case ROTATIONAL:
		{
			label face0 = getConsistentRotationFace(half0Ctrs);
			label face1 = getConsistentRotationFace(half1Ctrs);

			n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
			n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
			n0 /= mag(n0) + VSMALL;
			n1 /= mag(n1) + VSMALL;
			break;
		}
		default:
		{
			// Assumes that cyclic is planar. This is also the initial
			// condition for patches without faces.

			// Determine the face with max area on both halves. These
			// two faces are used to determine the transformation tensors
			label max0I = findMaxArea(pp.points(), half0Faces);
			n0 = half0Faces[max0I].normal(pp.points());
			n0 /= mag(n0) + VSMALL;

			label max1I = findMaxArea(pp.points(), half1Faces);
			n1 = half1Faces[max1I].normal(pp.points());
			n1 /= mag(n1) + VSMALL;
		}
	}

	if (mag(n0 & n1) < 1 - polyPatch::matchTol_())
	{
		if (debug)
		{
			Pout<< "cyclicPolyPatch::getCentresAndAnchors : Rotation :"
				<< " n0:" << n0 << " n1:" << n1 << endl;
		}

		// Rotation (around origin)
		const tensor reverseT(rotationTensor(n0, -n1));

		// Rotation
		forAll(half0Ctrs, faceI)
		{
			half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
			anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
		}
	}
	else
	{
		// Parallel translation. Get average of all used points.

		primitiveFacePatch half0(half0Faces, pp.points());
		const pointField& half0Pts = half0.localPoints();
		const point ctr0(sum(half0Pts)/half0Pts.size());

		primitiveFacePatch half1(half1Faces, pp.points());
		const pointField& half1Pts = half1.localPoints();
		const point ctr1(sum(half1Pts)/half1Pts.size());

		if (debug)
		{
			Pout<< "cyclicPolyPatch::getCentresAndAnchors : Translation :"
				<< " n0:" << n0 << " n1:" << n1
				<< " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
		}

		half0Ctrs += ctr1 - ctr0;
		anchors0 += ctr1 - ctr0;
	}

	// Calculate typical distance per face
	tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
}


// Calculates faceMap and rotation. Returns true if all ok.
bool Foam::cyclicPolyPatch::matchAnchors
(
	const bool report,
	const primitivePatch& pp,
	const labelList& half0ToPatch,
	const pointField& anchors0,

	const labelList& half1ToPatch,
	const faceList& half1Faces,
	const labelList& from1To0,

	const scalarField& tols,

	labelList& faceMap,
	labelList& rotation
) const
{
	// Set faceMap such that half0 faces get first and corresponding half1
	// faces last.

	forAll(half0ToPatch, half0FaceI)
	{
		// Label in original patch
		label patchFaceI = half0ToPatch[half0FaceI];

		faceMap[patchFaceI] = half0FaceI;

		// No rotation
		rotation[patchFaceI] = 0;
	}

	bool fullMatch = true;

	forAll(from1To0, half1FaceI)
	{
		label patchFaceI = half1ToPatch[half1FaceI];

		// This face has to match the corresponding one on half0.
		label half0FaceI = from1To0[half1FaceI];

		label newFaceI = half0FaceI + pp.size()/2;

		faceMap[patchFaceI] = newFaceI;

		// Rotate patchFaceI such that its f[0] aligns with that of
		// the corresponding face
		// (which after shuffling will be at position half0FaceI)

		const point& wantedAnchor = anchors0[half0FaceI];

		rotation[newFaceI] = getRotation
		(
			pp.points(),
			half1Faces[half1FaceI],
			wantedAnchor,
			tols[half1FaceI]
		);

		if (rotation[newFaceI] == -1)
		{
			fullMatch = false;

			if (report)
			{
				const face& f = half1Faces[half1FaceI];
				SeriousErrorIn
				(
					"cyclicPolyPatch::matchAnchors(..)"
				)   << "Patch " << name() << " : "
					<< "Cannot find point on face " << f
					<< " with vertices:"
					<< IndirectList<point>(pp.points(), f)()
					<< " that matches point " << wantedAnchor
					<< " when matching the halves of cyclic patch " << name()
					<< endl
					<< "Continuing with incorrect face ordering from now on!"
					<< endl;
			}
		}
	}

	return fullMatch;
}


Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace
(
	const pointField& faceCentres
) const
{
	const scalarField magRadSqr =
		magSqr((faceCentres - rotationCentre_) ^ rotationAxis_);
	scalarField axisLen = (faceCentres - rotationCentre_) & rotationAxis_;
	axisLen = axisLen - min(axisLen);
	const scalarField magLenSqr = magRadSqr + axisLen*axisLen;

	label rotFace = -1;
	scalar maxMagLenSqr = -GREAT;
	scalar maxMagRadSqr = -GREAT;
	forAll(faceCentres, i)
	{
		if (magLenSqr[i] >= maxMagLenSqr)
		{
			if (magRadSqr[i] > maxMagRadSqr)
			{
				rotFace = i;
				maxMagLenSqr = magLenSqr[i];
				maxMagRadSqr = magRadSqr[i];
			}
		}
	}

	if (debug)
	{
		Info<< "getConsistentRotationFace(const pointField&)" << nl
			<< "    rotFace = " << rotFace << nl
			<< "    point =  " << faceCentres[rotFace] << endl;
	}

	return rotFace;
}


// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //

Foam::cyclicPolyPatch::cyclicPolyPatch
(
	const word& name,
	const label size,
	const label start,
	const label index,
	const polyBoundaryMesh& bm
)
:
	coupledPolyPatch(name, size, start, index, bm),
	coupledPointsPtr_(nullptr),
	coupledEdgesPtr_(nullptr),
	featureCos_(0.9),
	transform_(UNKNOWN),
	rotationAxis_(vector::zero),
	rotationCentre_(point::zero),
	rotationAngle_(0),
	separationVector_(vector::zero)
{}


Foam::cyclicPolyPatch::cyclicPolyPatch
(
	const word& name,
	const dictionary& dict,
	const label index,
	const polyBoundaryMesh& bm
)
:
	coupledPolyPatch(name, dict, index, bm),
	coupledPointsPtr_(nullptr),
	coupledEdgesPtr_(nullptr),
	featureCos_(0.9),
	transform_(UNKNOWN),
	rotationAxis_(vector::zero),
	rotationCentre_(point::zero),
	rotationAngle_(0),
	separationVector_(vector::zero)
{
	dict.readIfPresent("featureCos", featureCos_);

	if (dict.found("transform"))
	{
		transform_ = transformTypeNames.read(dict.lookup("transform"));

		switch (transform_)
		{
			case ROTATIONAL:
			{
				dict.lookup("rotationAxis") >> rotationAxis_;
				dict.lookup("rotationCentre") >> rotationCentre_;
				rotationAngle_ = readScalar(dict.lookup("rotationAngle"));

				// Check rotation axis
				if (mag(rotationAxis_) < SMALL)
				{
					FatalErrorIn
					(
					    "cyclicPolyPatch::cyclicPolyPatch\n"
					    "(\n"
					    "    const word& name,\n"
					    "    const dictionary& dict,\n"
					    "    const label index,\n"
					    "    const polyBoundaryMesh& bm\n"
					    ")"
					)   << "Incorrect rotation axis: " << rotationAxis_
					    << abort(FatalError);
				}

				rotationAxis_ /= mag(rotationAxis_);

				// Check rotation axis if appropriate

				break;
			}
			case TRANSLATIONAL:
			{
				dict.lookup("separationVector") >> separationVector_;
				break;
			}
			default:
			{
				// no additional info required
			}
		}
	}
}


Foam::cyclicPolyPatch::cyclicPolyPatch
(
	const cyclicPolyPatch& pp,
	const polyBoundaryMesh& bm,
	const label index,
	const label newSize,
	const label newStart
)
:
	coupledPolyPatch(pp, bm, index, newSize, newStart),
	coupledPointsPtr_(nullptr),
	coupledEdgesPtr_(nullptr),
	featureCos_(pp.featureCos_),
	transform_(pp.transform_),
	rotationAxis_(pp.rotationAxis_),
	rotationCentre_(pp.rotationCentre_),
	rotationAngle_(pp.rotationAngle_),
	separationVector_(pp.separationVector_)
{}


Foam::cyclicPolyPatch::cyclicPolyPatch
(
	const cyclicPolyPatch& pp
)
:
	coupledPolyPatch(pp),
	coupledPointsPtr_(nullptr),
	coupledEdgesPtr_(nullptr),
	featureCos_(pp.featureCos_),
	transform_(pp.transform_),
	rotationAxis_(pp.rotationAxis_),
	rotationCentre_(pp.rotationCentre_),
	rotationAngle_(pp.rotationAngle_),
	separationVector_(pp.separationVector_)
{}


Foam::cyclicPolyPatch::cyclicPolyPatch
(
	const cyclicPolyPatch& pp,
	const polyBoundaryMesh& bm
)
:
	coupledPolyPatch(pp, bm),
	coupledPointsPtr_(nullptr),
	coupledEdgesPtr_(nullptr),
	featureCos_(pp.featureCos_),
	transform_(pp.transform_),
	rotationAxis_(pp.rotationAxis_),
	rotationCentre_(pp.rotationCentre_),
	rotationAngle_(pp.rotationAngle_),
	separationVector_(pp.separationVector_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::cyclicPolyPatch::~cyclicPolyPatch()
{
	deleteDemandDrivenData(coupledPointsPtr_);
	deleteDemandDrivenData(coupledEdgesPtr_);
}



// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::cyclicPolyPatch::initAddressing()
{
	polyPatch::initAddressing();
}


void Foam::cyclicPolyPatch::calcAddressing()
{
	polyPatch::calcAddressing();
}


void Foam::cyclicPolyPatch::initGeometry()
{
	polyPatch::initGeometry();
}


void Foam::cyclicPolyPatch::calcGeometry()
{
	polyPatch::calcGeometry();

	calcTransforms();
}


void Foam::cyclicPolyPatch::initMovePoints(const pointField& p)
{
	polyPatch::initMovePoints(p);
}


void Foam::cyclicPolyPatch::movePoints(const pointField& p)
{
	polyPatch::movePoints(p);
	calcTransforms();
}


void Foam::cyclicPolyPatch::initUpdateMesh()
{
	polyPatch::initUpdateMesh();
}


void Foam::cyclicPolyPatch::updateMesh()
{
	polyPatch::updateMesh();
	deleteDemandDrivenData(coupledPointsPtr_);
	deleteDemandDrivenData(coupledEdgesPtr_);
}


const Foam::edgeList& Foam::cyclicPolyPatch::coupledPoints() const
{
	if (!coupledPointsPtr_)
	{
		// Look at cyclic patch as two halves, A and B.
		// Now all we know is that relative face index in halfA is same
		// as coupled face in halfB and also that the 0th vertex
		// corresponds.

		// From halfA point to halfB or -1.
		labelList coupledPoint(nPoints(), -1);

		for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
		{
			const face& fA = localFaces()[patchFaceA];

			forAll(fA, indexA)
			{
				label patchPointA = fA[indexA];

				if (coupledPoint[patchPointA] == -1)
				{
					const face& fB = localFaces()[patchFaceA + size()/2];

					label indexB = (fB.size() - indexA) % fB.size();

					// Filter out points on wedge axis
					if (patchPointA != fB[indexB])
					{
					    coupledPoint[patchPointA] = fB[indexB];
					}
				}
			}
		}

		coupledPointsPtr_ = new edgeList(nPoints());
		edgeList& connected = *coupledPointsPtr_;

		// Extract coupled points.
		label connectedI = 0;

		forAll(coupledPoint, i)
		{
			if (coupledPoint[i] != -1)
			{
				connected[connectedI++] = edge(i, coupledPoint[i]);
			}
		}

		connected.setSize(connectedI);

		if (debug)
		{
			OFstream str
			(
				boundaryMesh().mesh().time().path()
			   /"coupledPoints.obj"
			);
			label vertI = 0;

			Pout<< "Writing file " << str.name() << " with coordinates of "
				<< "coupled points" << endl;

			forAll(connected, i)
			{
				const point& a = points()[meshPoints()[connected[i][0]]];
				const point& b = points()[meshPoints()[connected[i][1]]];

				str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl;
				str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl;
				vertI += 2;

				str<< "l " << vertI-1 << ' ' << vertI << nl;
			}
		}

		// Remove any addressing calculated for the coupled edges calculation
		const_cast<primitivePatch&>(static_cast<const primitivePatch&>(*this))
			.clearOut();
	}
	return *coupledPointsPtr_;
}


const Foam::edgeList& Foam::cyclicPolyPatch::coupledEdges() const
{
	if (!coupledEdgesPtr_)
	{
		// Build map from points on halfA to points on halfB.
		const edgeList& pointCouples = coupledPoints();

		Map<label> aToB(2*pointCouples.size());

		forAll(pointCouples, i)
		{
			const edge& e = pointCouples[i];

			aToB.insert(e[0], e[1]);
		}

		// Map from edge on half A to points (in halfB indices)
		EdgeMap<label> edgeMap(nEdges());

		for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
		{
			const labelList& fEdges = faceEdges()[patchFaceA];

			forAll(fEdges, i)
			{
				label edgeI = fEdges[i];

				const edge& e = edges()[edgeI];

				// Convert edge end points to corresponding points on halfB
				// side.
				Map<label>::const_iterator fnd0 = aToB.find(e[0]);
				if (fnd0 != aToB.end())
				{
					Map<label>::const_iterator fnd1 = aToB.find(e[1]);
					if (fnd1 != aToB.end())
					{
					    edgeMap.insert(edge(fnd0(), fnd1()), edgeI);
					}
				}
			}
		}

		coupledEdgesPtr_ = new edgeList(nEdges()/2);
		edgeList& coupledEdges = *coupledEdgesPtr_;
		label coupleI = 0;

		for (label patchFaceB = size()/2; patchFaceB < size(); patchFaceB++)
		{
			const labelList& fEdges = faceEdges()[patchFaceB];

			forAll(fEdges, i)
			{
				label edgeI = fEdges[i];

				const edge& e = edges()[edgeI];

				// Look up halfA edge from HashTable.
				EdgeMap<label>::iterator iter = edgeMap.find(e);

				if (iter != edgeMap.end())
				{
					label halfAEdgeI = iter();

					// Store correspondence. Filter out edges on wedge axis.
					if (halfAEdgeI != edgeI)
					{
					    coupledEdges[coupleI++] = edge(halfAEdgeI, edgeI);
					}

					// Remove so we build unique list
					edgeMap.erase(iter);
				}
			}
		}

		coupledEdges.setSize(coupleI);


		// Some checks

		forAll(coupledEdges, i)
		{
			const edge& e = coupledEdges[i];

			if (e[0] == e[1] || e[0] < 0 || e[1] < 0)
			{
				FatalErrorIn("cyclicPolyPatch::coupledEdges() const")
					<< "Problem : at position " << i
					<< " illegal couple:" << e
					<< abort(FatalError);
			}
		}

		if (debug)
		{
			OFstream str
			(
				boundaryMesh().mesh().time().path()
			   /"coupledEdges.obj"
			);
			label vertI = 0;

			Pout<< "Writing file " << str.name() << " with centres of "
				<< "coupled edges" << endl;

			forAll(coupledEdges, i)
			{
				const edge& e = coupledEdges[i];

				const point& a = edges()[e[0]].centre(localPoints());
				const point& b = edges()[e[1]].centre(localPoints());

				str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl;
				str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl;
				vertI += 2;

				str<< "l " << vertI-1 << ' ' << vertI << nl;
			}
		}

		// Remove any addressing calculated for the coupled edges calculation
		const_cast<primitivePatch&>(static_cast<const primitivePatch&>(*this))
			.clearOut();
	}
	return *coupledEdgesPtr_;
}


void Foam::cyclicPolyPatch::initOrder(const primitivePatch& pp) const
{}


// Return new ordering. Ordering is -faceMap: for every face index
// the new face -rotation:for every new face the clockwise shift
// of the original face. Return false if nothing changes (faceMap
// is identity, rotation is 0)
bool Foam::cyclicPolyPatch::order
(
	const primitivePatch& pp,
	labelList& faceMap,
	labelList& rotation
) const
{
	faceMap.setSize(pp.size());
	faceMap = -1;

	rotation.setSize(pp.size());
	rotation = 0;

	if (pp.empty())
	{
		// No faces, nothing to change.
		return false;
	}

	if (pp.size() & 1)
	{
		FatalErrorIn("cyclicPolyPatch::order(..)")
			<< "Size of cyclic " << name() << " should be a multiple of 2"
			<< ". It is " << pp.size() << abort(FatalError);
	}

	label halfSize = pp.size()/2;

	// Supplied primitivePatch already with new points.
	// Cyclics are limited to one transformation tensor
	// currently anyway (i.e. straight plane) so should not be too big a
	// problem.


	// Indices of faces on half0
	labelList half0ToPatch;
	// Indices of faces on half1
	labelList half1ToPatch;


	// 1. Test if already correctly oriented by starting from trivial ordering.
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	half0ToPatch = identity(halfSize);
	half1ToPatch = half0ToPatch + halfSize;

	// Get faces
	faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
	faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));

	// Get geometric quantities
	pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
	scalarField tols;
	getCentresAndAnchors
	(
		pp,
		half0Faces,
		half1Faces,
		ppPoints,
		half0Ctrs,
		half1Ctrs,
		anchors0,
		tols
	);

	// Geometric match of face centre vectors
	labelList from1To0;
	bool matchedAll = matchPoints
	(
		half1Ctrs,
		half0Ctrs,
		tols,
		false,
		from1To0
	);

	if (debug)
	{
		Pout<< "cyclicPolyPatch::order : test if already ordered:"
			<< matchedAll << endl;
	}


	// 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.)
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	if (!matchedAll)
	{
		label faceI = 0;
		for (label i = 0; i < halfSize; i++)
		{
			half0ToPatch[i] = faceI++;
			half1ToPatch[i] = faceI++;
		}

		// And redo all matching
		half0Faces = UIndirectList<face>(pp, half0ToPatch);
		half1Faces = UIndirectList<face>(pp, half1ToPatch);

		getCentresAndAnchors
		(
			pp,
			half0Faces,
			half1Faces,
			ppPoints,
			half0Ctrs,
			half1Ctrs,
			anchors0,
			tols
		);

		// Geometric match of face centre vectors
		matchedAll = matchPoints
		(
			half1Ctrs,
			half0Ctrs,
			tols,
			false,
			from1To0
		);

		if (debug)
		{
			Pout<< "cyclicPolyPatch::order : test if pairwise ordered:"
				<< matchedAll << endl;
		}
	}


	// 3. Baffles(coincident faces) converted into cyclics (e.g. jump)
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	if (!matchedAll)
	{
		label baffleI = 0;

		forAll(pp, faceI)
		{
			const face& f = pp.localFaces()[faceI];
			const labelList& pFaces = pp.pointFaces()[f[0]];

			label matchedFaceI = -1;

			forAll(pFaces, i)
			{
				label otherFaceI = pFaces[i];

				if (otherFaceI > faceI)
				{
					const face& otherF = pp.localFaces()[otherFaceI];

					// Note: might pick up two similar oriented faces
					//       (but that is illegal anyway)
					if (f == otherF)
					{
					    matchedFaceI = otherFaceI;
					    break;
					}
				}
			}

			if (matchedFaceI != -1)
			{
				half0ToPatch[baffleI] = faceI;
				half1ToPatch[baffleI] = matchedFaceI;
				baffleI++;
			}
		}

		if (baffleI == halfSize)
		{
			// And redo all matching
			half0Faces = UIndirectList<face>(pp, half0ToPatch);
			half1Faces = UIndirectList<face>(pp, half1ToPatch);

			getCentresAndAnchors
			(
				pp,
				half0Faces,
				half1Faces,
				ppPoints,
				half0Ctrs,
				half1Ctrs,
				anchors0,
				tols
			);

			// Geometric match of face centre vectors
			matchedAll = matchPoints
			(
				half1Ctrs,
				half0Ctrs,
				tols,
				false,
				from1To0
			);

			if (debug)
			{
				Pout<< "cyclicPolyPatch::order : test if baffles:"
					<< matchedAll << endl;
			}
		}
	}


	// 4. Automatic geometric ordering
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	if (!matchedAll)
	{
		// Split faces according to feature angle or topology
		bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);

		if (!okSplit)
		{
			// Did not split into two equal parts.
			return false;
		}

		// And redo all matching
		half0Faces = UIndirectList<face>(pp, half0ToPatch);
		half1Faces = UIndirectList<face>(pp, half1ToPatch);

		getCentresAndAnchors
		(
			pp,
			half0Faces,
			half1Faces,
			ppPoints,
			half0Ctrs,
			half1Ctrs,
			anchors0,
			tols
		);

		// Geometric match of face centre vectors
		matchedAll = matchPoints
		(
			half1Ctrs,
			half0Ctrs,
			tols,
			false,
			from1To0
		);

		if (debug)
		{
			Pout<< "cyclicPolyPatch::order : automatic ordering result:"
				<< matchedAll << endl;
		}
	}


	if (!matchedAll || debug)
	{
		// Dump halves
		fileName nm0(name()+"_half0_faces.obj");
		Pout<< "cyclicPolyPatch::order : Writing half0"
			<< " faces to OBJ file " << nm0 << endl;
		writeOBJ(nm0, half0Faces, pp.points());

		fileName nm1(name()+"_half1_faces.obj");
		Pout<< "cyclicPolyPatch::order : Writing half1"
			<< " faces to OBJ file " << nm1 << endl;
		writeOBJ(nm1, half1Faces, pp.points());

		OFstream ccStr
		(
			boundaryMesh().mesh().time().path()
		   /name() + "_faceCentres.obj"
		);
		Pout<< "cyclicPolyPatch::order : "
			<< "Dumping currently found cyclic match as lines between"
			<< " corresponding face centres to file " << ccStr.name()
			<< endl;

		// Recalculate untransformed face centres
		//pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
		label vertI = 0;

		forAll(half1Ctrs, i)
		{
			if (from1To0[i] != -1)
			{
				// Write edge between c1 and c0
				//const point& c0 = rawHalf0Ctrs[from1To0[i]];
				const point& c0 = half0Ctrs[from1To0[i]];
				const point& c1 = half1Ctrs[i];
				writeOBJ(ccStr, c0, c1, vertI);
			}
		}
	}


	if (!matchedAll)
	{
		SeriousErrorIn
		(
			"cyclicPolyPatch::order"
			"(const primitivePatch&, labelList&, labelList&) const"
		)   << "Patch " << name() << " : "
			<< "Cannot match vectors to faces on both sides of patch "
			<< name() << endl
			<< ".  Perhaps your faces do not match?"
			<< " The obj files written contain the current match." << endl
			<< "    Continuing with incorrect face ordering from now on!"
			<< endl;

			return false;
	}


	// Set faceMap such that half0 faces get first and corresponding half1
	// faces last.
	matchAnchors
	(
		true,                   // report if anchor matching error
		pp,
		half0ToPatch,
		anchors0,
		half1ToPatch,
		half1Faces,
		from1To0,
		tols,
		faceMap,
		rotation
	);

	// Return false if no change neccesary, true otherwise.

	forAll(faceMap, faceI)
	{
		if (faceMap[faceI] != faceI || rotation[faceI] != 0)
		{
			return true;
		}
	}

	return false;
}


void Foam::cyclicPolyPatch::syncOrder() const
{}


void Foam::cyclicPolyPatch::write(Ostream& os) const
{
	polyPatch::write(os);
	os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl;
	switch (transform_)
	{
		case ROTATIONAL:
		{
			os.writeKeyword("transform") << transformTypeNames[ROTATIONAL]
				<< token::END_STATEMENT << nl;
			os.writeKeyword("rotationAxis") << rotationAxis_
				<< token::END_STATEMENT << nl;
			os.writeKeyword("rotationCentre") << rotationCentre_
				<< token::END_STATEMENT << nl;
			os.writeKeyword("rotationAngle") << rotationAngle_
				<< token::END_STATEMENT << nl;
			break;
		}
		case TRANSLATIONAL:
		{
			os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL]
				<< token::END_STATEMENT << nl;
			os.writeKeyword("separationVector") << separationVector_
				<< token::END_STATEMENT << nl;
			break;
		}
		default:
		{
			// no additional info to write
		}
	}
}


// ************************************************************************* //
